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Abstract 



I We present a new computational methodology for the investigation 

. of gel electrophoresis of polyelectrolytes. We have developed the method 

Q ' initially to incorporate sliding motion of tight parts of a polymer pulled 

O , by an electric field into the bond fluctuation method (BFM). Such motion 

due to tensile force over distances much larger than the persistent length 
is realized by non-local movement of a slack monomer at an either end 
^ , of the tight part. The latter movement is introduced stochastically. This 

' new BFM overcomes the well-known difficulty in the conventional BFM 

^Nj ' that polymers are trapped by gel fibers in relatively large fields. At the 

same time it also reproduces properly equilibrium properties of a poly- 
mer in a vanishing filed limit. The new BFM thus turns out an efficient 
\ computational method to study gel electrophoresis in a wide range of the 

OP . electric field strength. 

J3 ■ 1 Introduction 

, Gel electrophoresis is a method which separates polyelectrolytes such as DNA 

' according to their length. This technique is an application of the phenomenon 

O , which polyelectrolytes exhibit different migration velocities in gel when they 

are pulled by an external electric field. Although in recent years many sophis- 
ticated methods such as pulsed-field techniques have been developed 0, ||, its 
underlying physics, even merely on steady field techniques, are not thoroughly 
understood. 

' The concept of reptation, which was proposed by de Gennes ||^, has been 

incorporated into theoretical analysis of gel electrophoresis. For example, the 
biased reptation model has succeeded to explain an empirical law of steady field 
experiments in the small field limit |5|, i.e., DNA mobility is proportional to 
reciprocal of its length However such predictions by the reptation theory 
are only on properties associated with dynamics of averaged conformations of 
DNA. 
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It is now well known that DNA dynamics in electrophoresis is more compli- 
cated and interesting. The molecular-dynamics (MD) simulation by Deutsch [Q 
on a freely jointed chain in a two-dimensional space with obstacles substituted 
for gel fibers first demonstrated that the chain migrates through obstacles tak- 
ing extended and collapsed conformations alternatively. Since this oscillatory 
behavior appears as a steady state the averaged conformation which the 
reptation model predicts is not stable. Also experimentally, since fluorescent 
microscopy has been invented and improved, evidences showing this instability 
have increased as well, and now details of DNA motion itself are of current 



interest in the study of gel electrophoresis ||9|, |10|, [11] y_2|. Among them the 
inch-worm like motion may be the most typical example which is observed for 
large DNA |l|. 

In studying polymer dynamics, numerical simulations have played important 
roles. One of them is the MD method, whose example was already mentioned 
above. Another powerful simulation is the bond fluctuation method (BFM), 
which is a Monte Carlo method utilizing description of a polymer(s) on a lattice 
[ p!^ . Its efficiency in examining gel electrophoresis in a small field was already 
demonstrated, but at the same time, its difficulty when applied to the case in 
a large field was pointed out The difhculty is that it takes huge MC steps 
(mcs's) for polymers once hanging on obstacles to get rid of them. As a result 
mobility in a large field becomes significantly smaller than in a relatively small 
field. This is an artifact of the numerical method since in actual experiments 
the mobility is observed to increase monotonically with increasing field ||l^ . A 
reason of this difficulty is considered that tensile force between monomers, which 
are of fundamental importance when each part of the chain is in nearly straight 
conformations, is not taken into account in the conventional BFM (c-BFM). 

The purpose of the present paper is to report a new BFM (n-BFM) which we 
have developed to overcome the above-mentioned difficulty by introducing to the 
c-BFM new stochastic processes which simulate sliding motions caused by tensile 
force. The new method turns out to be able to reproduce qualitative aspects 
of gel electrophoresis phenomena in a wide range of the field. In large fields 
polymer motion is free from trapping by obstacles and mobility monotonically 
increases (and tends to saturate) with increasing field as we have expected. 
In small fields, on the other hand, the results simulated by the new method 
coincide with those simulated by the c-BFM if the time unit is properly scaled. 
This means that, in this field range, the stochastic process newly introduced 
here contributes to the entropy effect due to polymer conformations similarly 
to what the c-BFM process does. Furthermore, in a certain limited range of 
field we have observed extended and contracted conformations of a polymer 
very frequently, though not quasi-periodically. The n-BFM is considered very 
efficient, complementary to the MD method, in studying gel electrophoresis. 

In the next section we describe the n-BFM in detail. One further detail in 
algorithm is explained in Appendix The results simulated by the n-BFM are 
presented and compared with those obtained by the c-BFM in Section ||, and 
the last section is devoted to concluding remarks of the present work. 
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2 Method 



In the present work we consider a, d = 2 square lattice version of the bond 
fluctuation method (BFM) to simulate gel electrophoresis. In this method 
each monomer is represented by a unit cell (4 lattice sites) and each bond has a 
variable length I but with the restriction 2 < I < ViS (see Fig. |l] in Appendix 
^). No monomer can come to any site which is already occupied by other 
monomers and by obstacles. The latter are substituted for gel fibers and each 
of them consists also of a unit cell. In the present work they are distributed 
periodically in both x and y directions with period of alattice units. A uniform, 
steady electric field is applied in the (1,1) direction. We use dimensionless 
electric field strength E = qSe/kgT, where q, £, e, kg and T denote charge of 
a monomer, bare electric field strength, lattice constant (~ persistent length, 
and is put unity in the present work), Boltzmann constant and temperature, 
respectively. The MC processes in the BFM consist of local random walks (of 
unit length) of each monomers in a potential due to the electric field. The 
excluded volume effect is satisfied and no bond crossing occurs in this BFM, 
which we call the conventional BFM (c-BFM). Actual simulations are done on 
square lattices of sizes 3M x 3M with periodic boundary conditions. Here M 
is the number of monomers in the chain. 

Although the c-BFM above described is shown its efficiency in examining gel 
electrophoresis in smaller fields, it has serious difficulty in larger fields p5|. The 
latter is easily understood from the inspection of a case shown in Fig. Ilia where 
we show a schematic picture of the U-shaped conformation of a chain hooked by 
an obstacle and pulled by a large downward field. Within the c-BFM the chain 
can escape from the trap only by the following process: a relatively slack part 
of the chain, created around the end segment, climbs up sequentially against a 
potential slope of the field up to the position of the obstacle. The probability for 
this process to occur is the smaller for the larger field. For example, the mobility 
/i of the M = 200 chain simulated by the c-BFM increases with increasing E 
up to ~ 0.01, but it starts to decrease and tends to vanish when E is further 
increased. Here the mobility fi is simply evaluated by 

^^(XoM^XoM/^, (1) 

tf - u 

in our units, where (• • ■ ) indicates the average over the MC runs. In Eq. (|l|) Xq 
is the displacement of center of mass in the field direction. The starting time 
ti of observation is chosen to eliminate influence of an initial conformation, and 
tf is chosen by the condition that {XQ{tf) — XQ{ti)) is larger than c^a or c^i?i 
with > 10, where Ri is radius of gyration of the chain. 

In actual experiments DNA can slide off an obstacle even if it is temporarily 
trapped by the latter. For example, Volkmuth and Austin fabricated quasi-two- 
dimensional microlithographic array of posts and looked at a sequence of moving 
images of 100 kb DNA in steady field of l.Ov/cm (corresponding to i? ~ 0.005) 
[ p!7| . They observed that the DNA hooks one of the posts, is extended to nearly 
its full contour length, and then slides from the post. This indicates that tensile 
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force plays an important role in such sliding process because the chain is fully 
extended in the process. 

It is then natural to introduce to the c-BFM the following non-local processes 
which simulate realistic sliding motions of DNA due to tensile force. In case of 
the U-shaped conformation of Fig. |l|a, we simply remove the end monomer of 
the shorter arm and join it to the other end of the chain. For the M-shaped 
conformation shown in Fig. |l|b, we remove one of the monomers in a dangling 
part at the center to the end of the longer arm. 

Starting from the above intuitive idea, we have constructed an algorithm 
in which non-local dynamics of monomers is systematically and stochastically 
introduced. For this purpose we first define monomers, on which such non-local 
processes are tried. They should be in a slack part of the chain, through which 
tensile force cannot transmit. We call them slack monomers (s- monomers). 
The s-monomers defined consist of monomers at both ends of a chain and of 
monomers whose nearest neighboring monomers are not separated from each 
others by distance larger than or equal to 4. Parts of the chain between neigh- 
boring s-monomers are regarded as tight parts, through which tensile force trans- 
mits. They are considered to slide either partially or as a whole depending on 
non-local movement of s-monomers which we next introduce. 

Our guiding principle to specify a procedure of non-local movement of s- 
monomers is to choose such stochastic processes that any s-monomer moves so 
as to fulfill detailed balance condition in equilibrium. This restriction, however, 
does not specify a procedure uniquely. Among possible procedures we adopt 
the following one: 

i) Choose one s-monomer [monomer in Figs. |^a and b shown for examples] . 

ii) Count the number of pairs n, on which we can try to move the chosen s- 
monomer, and which are in a region between its neighboring s-monomers 
including ends of the chain [pairs 12, 23, 34, and 45 for case 2a ('end' for 
case 2b) between monomers and 5 ('end'), and 1'2', 2'3', 3'4', and 4'5' 
between monomers and 5', and so n = 8]. 

iii) Choose randomly one of these pairs [34 ('end')] with the probability 1/n. 

iv) Choose randomly one of the allowed conformations putting a monomer 
between the chosen pair [3a4 (4a)] with the probability 1/W. By con- 
struction the moved monomer is an s-monomer. 

v) Accept the movement of the s-monomer [from s-monomer to s-monomer 
a] according to the weight w{/S.X) defined by 

^(^^) = e-^Ax+e^Ax ^ (2) 

where AAT is the displacement of the s-monomer in the (positive) direction 
of the field. 

In (iv) above, the value of W is chosen to be greater than max{Wi, .., W„} 
where Wi is the number of the allowed conformations putting a monomer to the 
i-th pair. In the present work we fix ly = 23 which is the possible maximum 
value of Wi when the i-th pair is one of the 'ends'. Then the probability of the 
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movement of the s- monomer to one specified conformation is given by 1/nW 
multiplied weight w{AX) , while that of its reversed process is 1/nW multiplied 
by weight ■w{—AX). This guarantees the detailed balance condition. It is noted 
here that the acceptance ratio of one non-local movement without specifying 
the moved conformation at all, rmovo, is given by rmove — Wi/2W in the limit 
EAX 0, where Wi is the average of Wi's including the further restriction 
below described. 

There remain, however, two problems which require further considerations. 
One is concerned with the detailed balance condition. By the procedure de- 
scribed above it may happen that the number of s-monomers changes depending 
on positions of the monomers around the chosen pair [in case Fig. ||a, besides 
the a-monomer, monomers 3 or 4 may become an s-monomer]. If this happens, 
the number n for the moved s-monomer also changes after the movement (ex- 
cept for the case that the movement involves an 'end' s-monomer). Then, by the 
above procedure the probability of the reverse process violates the detailed bal- 
ance condition. We get rid of this problem simply by accepting only processes, 
by which the number n for the moved s-monomer is conserved. 

The other problem remained is how to exclude possible bond crossing which 
may occur when an s-monomer is moved to an end of the chain. We are faced 
to the same problem when we try to construct an initial conformation of a 
chain by the self-avoiding random walk. We have developed an algorithm which 
recognize whether bond crossing occurs or not by referring to only local data. 
It is explained in Appendix |^. 

The new BFM (n-BFM) we propose is a combined process of non-local move- 
ments of s-monomers by the above procedure and the c-BFM process]] : each 
odd (even) MC steps of the c-BFM are followed by trials of non-local movements 
of odd (even) numbered s-monomers. 

3 Results 

Let us begin with comparison of equilibrium properties in i? = simulated by 
both the c-BFM and the n-BFM. As for a representative of static equilibrium 
properties we show four sets of data of radius of gyrations Ri in Fig. ^ in cases 
with obstacles (a — 20) and without obstacles (a = oo) obtained by the two 
BFM's. They all coincide with each others and fit well with the power law 
i?i cx with 1/ = 0.75 ± 0.02 which is consistent with the previous results 
[ p^ . There is no significant difference between the cases a = 20 and a = cx) in 
Ri-M dependences. This is because the concentration of obstacles is too small 
to change the value of exponent v. |2| . 

In Fig. ^ the diffusion constants Dq are plotted against M to compare 
behaviors of fluctuation simulated by the two BFM's. Here Dq is evaluated 

^For a local movement of monomers in the c-BFM in the present work, we have adopted 
the same weight ui(AX) with AX = itl/\/2 after choosing one of the four neighboring sites 
to move with equal probability. 
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simply by 



^^^2(^7^ with (Ai?G)' = ([RG(i/)-RG(tOf>, (2) 

where Rg is coordinate of the center of mass of the chain. The data actually 
used are those in the range Ri < Ai?G < c^Ri where cd is at least larger than 
2. In case without obstacles both BFM's yield Dq (x as expected, while 

in case a = 20 the same power law dependence Dq^° with v-q ~ 1.71 ± 0.02 
(n-BFM) and 1.73 ± 0.04 (c-BFM) is obtained. However Dq by the n-BFM is 
larger by factor 2^3 than that by the c-BFM. 

A more stringent check of the new algorithm may be whether it reproduces 
the fluctuation-dissipation theory, or the Einstein relation /z/M = Dq, where 
fi is the mobility in the limit E Q. We plot /i/Af for various E and Dq 
against M in Fig. |^. Clearly we see the Einstein relation holds well. The 
results described so far indicate that our new algorithm reproduces equilibrium 
properties of a polymer quite satisfactorily. 

Let us next examine stationary properties in finite E simulated. In order 
to compare those obtained by the two BFM's, we show the ratio /icon/M for 
various M with a = 20 plotted against E in Fig. ^. Here /Xcon and fi are 
mobilities simulated by the c-BFM and the n-BFM, respectively. In small fields 
E < 5.0 X 10~^ the ratio is almost constant and is also independent of M. The 
constant value is equal to the ratio of the diffusion constants Dq of the two 
BFM's mentioned above. 

The ratio /icon//^ deviates from the constant value when E exceeds a certain 
value Etii which depends on M. Its deviation is more significant for chains with 
larger M, and the ratio tends to vanish at E much larger than Eth- This result 
simply reflects the difficulty of the c-BFM described in the previous section and 
is considered to have nothing to do with the n-BFM. Actually there seems no 
crossover-like behavior around Eth in obtained by the n-BFM as shown in 
Fig. 0. One sees in the figure that /i of each M increases monotonically with 
increasing E, and tends to saturate at largest E (~ 0.1) we have simulated. Thus 
the data in Figs. ^ and |^ tell us that we have in fact overcome the difficulty of 
the c-BFM by means of the n-BFM we have proposed. 

In order to get further insights of the n-BFM which gives rise to quite sat- 
isfactory results so far demonstrated, we have examined some details of the 
s-monomer movement in a system with M = 100 and a = 20. Interestingly, the 
acceptance ratio of the s-monomer movement, rmovc, defined in the previous 
section is rather small (~ 3%) and does not depend sensitively on E in the 
range examined as shown in the inset of Fig. |^. On the other hand, the average 
number of the s-monomers, Ms, is rather large: Mg ~ 28 at S = and Ms — 24 
at E ^ 0.05. 

These results are interpreted as follows. In the field range investigated here 
{E < 0.05), the entropy effect dominates to determine conformations of a chain. 
The latter is therefore in a coiled conformation on average. The s-monomer 
movement in this situation is regarded as local fiuctuation of the conformation 
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due to the entropy effect rather than due to the tensile force. Although Tmovc is 
rather small, the s-monomer movements contribute to displacement of the center 
of mass of the chain which is comparable in magnitude with that by one mcs of 
the c-BFM process. This explains why Dq and /i by the n-BFM are larger than 
the corresponding quantities by the c-BFM (note that one mcs in the n-BFM 
consists of one mcs of the c-BFM process and trials of non-local movement on 
half of the s-monomers). Such efficiency of the s-monomer movement to Dq 
and fj, may be attributed to the fact that the s-monomer moves along the chain 
direction which is hardly realized by the c-BFM process. 

The main part of Fig. || is the histogram of the relative ratio of the s- 
monomer movements accepted, rAm, against the distance (along the chain). 
Am, they moved. Notice that the movement with Am ^ M/2 occurs even at 
E = 0, though TAm is quite small (the abscissa is in the logarithmic scale). This 
rAm &t E = reflects how frequently, in fact quite rarely, extended conforma- 
tions occur by fluctuation in the coiled conformation of the chain. 

An important observation of Fig. ^ is that rAm for E > 0.02 is much 
enhanced at large Am, though the absolute magnitude is still quite small. This 
means that even in this range of E extended conformations occur with low 
probability. But once they are present, the s-monomer movements with large 
Am occur with relatively high probability. This situation is just what we have 
first intuitively expected: the tensile force is considered to play an important 
role on dynamics of the chain. The low probability of such non-local movements 
of the s-monomer, which nevertheless overcomes the difficulty in the c-BFM 
and reproduces smooth, monotonic ii^-dependence of the mobility, justifies our 
introduction of the s-monomer movements as a stochastic process. 

In Fig. ^ we show the n versus log M plot using the same data shown in Fig. 
. It can be seen that /x for each E decreases monotonically as M increases, and 
that they are strongly dependent on M when E is small enough. As E increases, 
however, it gradually looses M-dependence from larger values of M, and finally 
11 becomes almost independent of M when E is quite strong {E ~ 0.1). These 
behavior of n is qualitatively consistent with experimental results [p^ . 

Finally we note that chain conformation dynamics simulated by the present 
n-BFM is actually complicated. Particularly in a certain limited range of 
E (~ 0.01) sequences of contraction and extension are observed. One of such 
sequences are shown in Fig. |l^ which are snapshots of conformation of an 
M — 200 chain in every 4 x lO** mcs in a certain MC run. A chain moving 
in a relatively collapsed form (a) is trapped by a certain obstacle (b), deforms 
to a V- or U-shaped conformation (c,d), slides off the obstacle (e), and then 
tends to form a collapsed conformation again (f). In contrast to experimental 
observation [ p^ , however, quasi-periodical alternation between contracted and 
extended conformations has not yet been observed, or periods between the two 
conformations are rather random. 
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4 Concluding Remarks 



By introducing stochastic, non-local movements of slack parts of the polymer 
(s- monomers) into the conventional bond fluctuation method (c-BFM) , we have 
constructed a new BFM (n-BFM) algorithm which overcomes the difficulty of 
the c-BFM applied to gel electrophoresis in relatively large field, namely, poly- 
mers once hooked by gel fibers become unable to get rid of them. In smaller 
fields, on the other hand, the new stochastic process gives rise to conformational 
change of the polymer which is interpreted as the ordinary entropy effect. The 
present n-BFM thus turns out to be able to reproduce qualitative aspects of gel 
electrophoresis phenomena in a wide range of the field. 

The n-BFM is considered more effective for a denser gel (a smaller a). Actu- 
ally the preliminary analysis of the mobility in case a = 12 yields /icon/A* — 0.25 
at small E, which is smaller than that of the case a = 20 shown in Fig. ^. This 
tells that non-local s-monomer movements make it easier for parts of polymer 
to escape obstacles, and so even the entropic conformation change in smaller 
fields is fastened as compared with the c-BFM. 

There are many problems to be explored further. Among them are, an 
improvement of algorithms to increase the acceptance ratio of the s-monomer 
movements, rmovc, without violating the detailed balance, and quantitative com- 
parisons with experimental observations. For the latter purpose we have to 
adopt model systems with more realistic distribution and/or shape of obstacles, 
and to extend the algorithm to ci = 3 system. These problems are now under 
investigations. 
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A Rejection procedure for bond crossing 

Here we explain a method how to reject a non-local movement of an s-monomer 
to one end of a chain which associates with bond crossing. For this purpose we 
use, as a simple example, an M = 7 BFM chain drawn in Fig. |T^. In addition to 
monomers and connecting bonds, we introduce a local function represented 
by small squares on each site n. This function contains information on position 
of each site n relative to nearby bonds which are now specified as vectors i, z -I- i. 
When takes a value 'lA' ('iB'), it means that site n is near the bond + i 
(both distances from monomers i and i -I- 1 are less than 4) and is in its wright 
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(left) hand side. Since site n can be near other bonds $(n) is, in general, 
multi-valued. If, on the other hand, site n is far from any bonds $(n) ='null'. 

By making use of the local function $(n) thus defined we can check bond 
crossing as follows. To be simple, let us consider the case removing end monomer 
1 and putting it to the other end monomer 7, as shown in Fig. |ll|. If monomer 
7 touches a site with $ —'2A\ the moved monomer should not touch sites with 
$ =2B\ Monomer a is such a case where bonds 27§ and 7, a cross and so 
is rejected. Other cases such as monomer (5 are all accepted from the present 
criterion on bond crossing. 
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Figure 1: Schematic illustration of (a) the U-shaped and (b) M-shaped confor- 
mations pulled down by the field. 




Figure 2: Examples of s-monomers and their movements, (a) The case in which 
no end monomer is involved, (b) The case in which an end monomer is involved. 
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Figure 3: Dependence of radius of gyration Rj on length M. 
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Figure 4: Dq-M plots for a = 20, oo simulated by the two BFMs. ud denotes 
the exponent of the power law dependences. 



11 



10-1 



10" 



Dg or ^ 



M 



10" 



10" 



10 



® 



Dg : B - 
IJ/M,E=2.0xlO~' --e 
£=5.0x10"* ^ A 
£=10x10"* ^v- 



£=20x10"* O 

i , 



20 30 40 50 



M 



100 



200 300 400 



Figure 5: Diffusion constant Dq and i^/M plotted against M. 




gurc 6: The mobility ratio yUcon/M plotted against E. Where /icon is obtained 
the c-BFM and /x by the n-BFM. 
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Figure 7: Plots of vs. log-E. 
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Figure 8: Histograms of the s-monomer movements accepted, rAm, against the 
distance (along the chain) Am for the chain with M = 100 and a = 20. In the 
inset is shown the acceptance ratio of the s-monomer movement, rmove- 
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Figure 9: Plots of log M against fi. 




Figure 10: Time development, (a), (b), . . . , (f), of M = 200 chain conformations 
in field E || (1, 1) and E = 1.02 x IQ-^. The snapshots at every 4.0 x lO^mcs 
are shown. Grid points represent obstacles (o = 30). 
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Figure 11: An example of conformations of a M = 7 BFM chain on a square 
lattice. The unit square with thick broken hues represents an obstacle. Small 
squares on sites denote the local function $(n). 
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